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ABSTRACT 

We have performed hydrodynamic simulations of galaxy formation in a cold dark matter 
(ACDM) universe. We have followed galaxy formation in a dark matter halo, chosen to have 
a relatively quiet recent merger history, using different models for star formation and feed¬ 
back. In all cases, we have adopted a multi-phase description of the interstellar medium and 
modelled star formation in quiescent and burst modes. We have explored two triggers for star- 
bursts: strong shocks and high gas density, allowing for the possibility that stars in the burst 
may form with a top-heavy initial mass function. We find that the final morphology of the 
galaxy is extremely sensitive to the modelling of star formation and feedback. Starting from 
identical initial conditions, galaxies spanning the entire range of Hubble types, with /i-hand 
disc-to-total luminosity ratios ranging from 0.2 to 0.9, can form in the same dark matter halo. 
Models in which starbursts are induced by high gas density (qualitatively similar to models in 
which feedback is produced by AGN) generate energetic winds and result in galaxies with an 
early-type morphology. Models in which the starbursts are induced by strong shocks lead to 
extended discs. In this case, the feedback associated with the bursts suppresses the collapse 
of baryons in small haloes, helping to create a reservoir of hot gas that is available for cooling 
after z ~ 1, following the bulk of the dynamical activity that builds up the halo. This gas then 
cools to form an extended, young stellar disc. 
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1 INTRODUCTION 

Understanding galaxy formation is a challenging problem whose 
solution will require a concerted approach combining observa¬ 
tional, semi-analytic and numerical work. There have been substan¬ 
tial advances on all these fronts in the past decade, but fundamental 
questions remain unanswered. One of the most troublesome is the 
inability to produce realistic spiral galaxies in numerical simula¬ 
tions that start from initial conditions appropriate to the concor¬ 
dance cold dark matter cosmology (ACDM) and assume prescrip¬ 
tions for gas and stellar processes analogous to those included in 
semi-analytic models, or inferred directly from observations. 

Pioneering N-body/gasdynamic al simu lations b y 
iKatz & Gun J <199ll) . iNavarro & Ben3 <199 ill and iKatzl < 1 9921) . 
already included many of the processes, which are generally 
regarded as essential for galaxy formation: radiative cooling, star 
formation, and feedback effects generated by energy released from 
associate d with supemova e (SNe) and stellar winds. The simula¬ 
tions b y |Navarro_&Bend (l99l|), a nd late r by INavarro & White! 
<1994 and INavarro, Frenk & White! < 1 99.4 . assumed CDM initial 
conditions. They uncovered a fundamental problem that prevents 
the formation of disc-dominated galaxies: the so-called “angular 
momentum problem.” This is a generic problem that can be 
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traced back to very nature of the hierarchical clustering process 
characteristic of s tructure growth from CDM initial conditions 
<Frenk et al Jj 198.4 . At early times, small, dense cold dark matter 
haloes form. Radiative cooling is very efficient and a large fraction 
of their gas, some of which will turn into stars, cools into their 
centres. As haloes merge to form larger objects, their orbital 
angular momentum is drained by dynamical friction and exported 
to the dark matter at the outskirts of the new haloes. Much of the 
original angular momentum of the baryonic material is lost and the 
resulting galaxies become too centrally concentrated. 

[Weil, Eke & Efstathioul 1 1998 ) and lEke. Efstathiou & Wrighl 
<2000) " showed that if cooling is artificially suppressed un¬ 
til the host haloes are well established, then the simula¬ 
tions can produce galaxies that are less centrally concentrated 
and have higher specific angular momenta. Two ways to pre¬ 
vent the early collapse of protogalactic clouds have been pro¬ 
posed: stronger feedback than that provided by_standard treat- 
ments o f supemovae JSommer- Larsen. Gelato. & Vedell 1 1 999< 
iThacker & CouchmarfeOOll) and a modification of the cosmologi¬ 
cal framework, replacing CDM with warm dark matter (which does 
not in duc e small-scale fluctuati ons: ISommer-Larsen & Dolsoxl 
<200lh and lGovernato et alJ<2004 ). There are also indications that 
spurious numerical effects, arising from the very nature of the 
smoothed particle hydrodynamics (SPH) technique used in the sim¬ 
ulations, cause transfer of angular momentum from the cold gas 
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disc to the hot halo gas, s ig nific a ntly con tributing to the angular 
momentum problem i lOkamoto et alfcoQ^l . Aside from numerical 
effects, it is clear that the angular momentum problem in simula¬ 
tions is telling us something quite fundamental about the nature of 
star formation and feedback processes in galaxy formation. 

Some recent simulations have yielded more 
promising disc galaxies in the CDM framework. 
ISommer-Larsen. Gotz & Portinaril ( 2002 '! were able to gener¬ 
ate a variety of morphological types, including discs, by assuming 
a self-propagating star formation model combined with very 
efficient SNe feedback. In their model, star formation proceeds 
in two modes, “early time” and “late time”. In the early time 
mode, star formation is very efficient and the feedback is very 
strong; in the late time mode, the star formation efficiency is low 
and there is no feedback. iGovernato et alJ 11200 4 > produced a disc 
galaxy employing standard star formation and feedback recipes 
and claimed that numerical resolution is the primary cause of the 
angular momentum problem. Abadi et al. (2003a, 2003b) were 
able to generate a galaxy resembling early-type spirals for which 
they calculated detailed photometric and dynamical properties. 
iRobertson et alJ (20041) adopted the multi- phase model for the 
star-fo rming interstellar medium (ISM) of lspringel‘& Hernquist 
(2003), which stabilises gaseous discs against Toomre’s instability, 
and produced a galaxy having an exponential surface brightness 
profile. Note that earlier work by ISteinmetz & Navarrcl (1999) 
follow ed galaxy fo rmati on in randomly chose n haloes, while 
lAbadi et*all l2003al) and IGovernato et al J (20041) selected haloes 
that have no major mergers at low redshift (z < 2 in Governato et 
al.). In the ACDM model, haloes with such quiet merger histories, 
while favourable for disc formation, have a number density today 
which is too low to account for the observed number density of 
spiral galaxies. 

While early simulations already inclu ded sim plified 
treatments of chemical evolution ( Steininetz & Mulled 119951 ; 
iRaiteri. Villata & Navarrolll996l : lBerczikl^9^) . more recent sim¬ 
ulations have considered Type la supemovae (SNe la) in addition 
to Type II supernovae (SNe II), sometimes relaxing the instanta¬ 
neous recycling approximation HRA) (Li a. Portinar i. & Carrara 
20021: iKawata & Gibsorj 120031: iKobavashil l2004) .' In particular, 
Kawata & Gibson ( 2003 ) have argued that the non-instantaneous 
nature of SNe II and SNe la isJmportant for the dynamical evolu¬ 
tion of galaxies. [Brook et~ . ( 2004 1 have shown that the chemical 
abundance of the halo stars could be an additional constraint for 
the modelling of feedback. 

In this paper, we present a new series of simulations of 
galaxy formation. The main difference with previous work is 
that we consider an unconventional star formation model which, 
however, has been claimed to be required to explain the prop¬ 
erties of high redshift sub-millimeter and Lyman-break galaxies 
(iBaugh et~alJ 200^)^ as well a s the_metallicity of the intraclus- 
ter medium NagashimaetalJ 12005 ah and of elliptical galaxies 
(Nagashima et al . 2005bh . In this model, star formation normally 
proceeds in a quiescent mode with a standard initial mass func¬ 
tion (IMF). However, when a major merger occurs, it triggers a 
burst of star formation with a top-heavy IMF. The distinction be¬ 
tween quiescent and burst modes is not particularl y cont roversial 
but the adoption of a top-heavy IMF is. lBaugh et alJ 12005 ) claimed 
that the number density of sub-millimeter galaxies, in particular, 
is i mposs ible to explain without a top-heavy IMF in bursts, while 
iNagashima et alJ (2005a) claimed that the ratio of alpha to iron 
peak elements in the intracluster medium can only be understood 
if bursts have a top-heavy IMF. A review of the observational and 


theoretical evidence for and against a top-heavy IMF may be found 
in these papers and references therein. 

Our main result is that the adoption of a top-heavy IMF in 
bursts eases the formation of a large disc primarily because at early 
times, when mergers are more important, more energy per unit of 
mass turned into stars is made available for feedback. However, we 
also find that varying the criteria for a burst or varying the IMF 
in bursts can have a strong effect on the final morphology of the 
simulated galaxy. 

This paper is organised as follows. We describe our simula¬ 
tion code and model for star formation and feedback in Section 2. 
We present the details of our simulations in Section 3 and the re¬ 
sults in Section 4. Finally, in Section 5, we present a discussion and 
summary of our conclusions. 


2 SIMULATION CODE 

We use the parallel PM-TreeSPH code GADGET2 (Springe! 
[2005|), a successor of the TreeSPH code GADGET 
(Snringel. Yoshida & White] EoOlt). The hy drody namics are 
solved using an SPH algorithm lLuc\ll977l : lGingold & Mra taghan 
Il977t) and a “conservative entropy” formulatio n that manifestly 
conserves energy and entropy iSpringel & Hemauisti2002l) . 

Here, we briefly describe modifications we have made to the 
code so as to adapt it for our purposes. In the conservative entropy 
formulation, the smoothing length of an SPH particle should be 
given by 



pi = M ngb = const., 


( 1 ) 


where hi and pi are the smoothing length and density of the ith 
SPH particle respectively, and M ng b represents the mass in its 
smoothing volume. We set M ng b = 40m O ri g throughout this pa¬ 
per, where m or i g is the original SPH particle mass. The default 
implementation of GADGET2 keeps the number of neighbouring 
particles in the smoothing volume constant, whereas we chose to 
keep the mass resolution constant as in equation Q. Since in our 
simulations the gas particle masses vary due to star formation and 
feedback, these two choices are not equivalent. We also adopt the 
phase decoupling technique introduced bv lOkamoto et alJ (; 2003t) to 
avoid spurious angular momentum transfer from cold gas discs to 
surrounding hot halo gas. We will describe the modelling of cool¬ 
ing, star formation, and feedback in the following subsections. 


2.1 Gas cooling 

We calculate the cooling/heating rate and ionization state of 
each particle by assuming collisional ionization equilibrium 
and the presence of an evolving but uniform UV background 
(llaardt & Madaii 11996 1 that is switched on at z = 6. Inverse 
Compton cooling is also included. In order to take the metallic- 
ity depe ndence in to account, w e use t he appropriate cooling tables 
given bv ISutherland & Dopital»199 3) at T > 10 4 K. The coolest 
gas in overdense regions typically has a temperature T ~ 10 4 K, 
because we do not include molecular cooling or metal cooling be¬ 
low 10 4 K. 

Metals are carried by particles and once assigned to a particle, 
they remain with it. Nonetheless, effective mixing takes place be¬ 
cause we use the smoothed metallicity (smoothed in the same way 
as other SPH quantities) when computing the cooling rates. 
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2.2 Star formation and feedback 

In general, star formation in numerical simulations is modelled as 

£ = ( 2 ) 

CLZ £dyn 

where p*, p g , tdyn = (47rGp g ) -1 ^ 2 , and c* are the stellar density, 
the gas density, the local dynamical time, and a dimensionless star 
formation efficiency parameter respectively. 


Typically, c« ~ 1/30 has been used in simulations of disc 


formation (e.g. Steinmetz & Navarro 1999 

Thacker & Couchman 

2001 

; Abadi et al. 2003a; Governato et al 

2004; Robertson et al 

2004 

) in order to reproduce the observed gas mass fraction and/or 


the Kennicutt (1998) law. Meza et al. (2003 

) showed that the same 

star formation model used by Abadi et al. 

2003aj) could also lead 


to the formation of an elliptical galaxy if the initial conditions 
were such that a major merger took place at low redshift. On the 
other hand, some simulations of elliptical galaxy formatio n, have 
used jarge values of c* = 0.5 ~ 1 te.g. tKawata & Gibsonl2004l : 
lKohavashl200i) . often with strong feedback, in order to reproduce 
the observed sizes and/or colours of elliptical galaxies. Sommer- 
Larsen et al. (2003) combined these two star formation recipes by 
introducing ‘early’ and Tate’ star formation modes. In the early star 
formation mode, they adopted c* = 1, together with a low thresh¬ 
old density for star formation and very strong feedback while, in the 
late star formation mode, they adopted c* = 1/40, and prevented 
any feedback. With this prescription, they were able to reproduce a 
wide range of morphological types in their simulations. 

We consider two physically motivated star formation modes, 
which we call ‘quiescent’ and ‘burst’. We model the quiescent 
mode as self-regulated star formation using a multi-phase ISM 
model based on that developed bv lspringel & Hernauisa (EoP-lh . but 
modified to avoid the instantaneous recycling approximation (IRA) 
in which short-lived stars assumed to die immediately as they form. 
We use the same model for the burst mode, except that we adopt a 
shorter star formation timescale (c* =0.5) and a flatter IMF with 
the intention of making self-regulation an unstable process. The 
motivation for these choices was discussed in § 1. 


2.2.1 A multi-phase model for star forming gas 

Following lSpringel & Hemquis t (2003), the dense ISM gas is pic¬ 
tured as a two-phase fluid consisting of cold clouds and an ambient 
hot phase whose energy is supplied by SNe explosions. We now 
briefly explain our model, including the modifications we have im¬ 
plemented. 

In[Snringel & Hernquistj 1200 3l) . the densities of cold clouds 
(p c ) and the hot phase (ph) are related to the total gas density (p) as 
p = p c -f- p h . By considering the volume that each phase occupies, 
we impose the following relations: 


Msph 

M c M h 

(3) 

P 

pc ph 

Msph = 

M c + Mh, 

(4) 


where Msph, M c , and Mh are the particle mass, mass in cold 
clouds, and mass in the hot phase associated with the particle. We 
further assume the cold clouds are in pressure equilibrium with the 
hot ambient phase, namely 

PMeff ~ PcUc — PhMh, (5) 

where u e g, is the effective specific internal energy of an SPH parti¬ 


cle, and it c and Mh are the specific internal energies of cold clouds 
and hot phase. 

We describe our multi-phase model, first assuming the IRA, 
in which star formation, cloud formation by thermal instability, 
evaporation of clouds, and heating of the hot phase by SN ex¬ 
plosions are inclu ded. This model is almost identical to that of 
iSpringel & Hemquisll 1200 3) except for the definition of the den¬ 
sity of each phase and the metallicity dependence. We will then 
show how to relax the IRA. 

Star formation takes place on a star formation timescale i*: 
dM* _ M c 

J, , 7 (6) 

at £* 

where M» is the stellar mass. The star formation timescale, /*, is 
taken to be compatible with equation J3 - for p > pth we have: 

/ \ -i/2 

f.(p) = *2{ —) , (7) 

\ Pth ) 

where the value of /* is chosen to match the Kennicutt law 
jKennicutill998! ) and pth is a threshold density above which SPH 
particles become multi-phase and are eligible to form stars. For 
P < Pth, no star formation takes place. 

The short-lived stars immediately die and return mass and re¬ 
lease energy, esn, for each solar mass in stars formed, as SNe II. 
The heating rate arising from SNe is therefore 


di 


(. MhUh ) 


: £SN- 


dtf, 
d t 


( 8 ) 


We assume that SN explosions also evaporate the cold clouds, 
and transfer gas back into the ambient hot phase. 


dM c 


d t 


- \ £SN Me 

BV MSN t * 


(9) 


Here, A is the efficiency parameter and the specific supernova en¬ 
ergy, msn = csn JJ 3 , where /3 is the mas s fraction of short-lived 
stars. Following tMcKee & Ostrike3 il977l ). the efficiency parame¬ 
ter has the density dependence 


A(p) = Aq 



( 10 ) 


Finally, we assume the cold clouds form and grow through 
thermal instability, that is 


dM c 


di 


dM h _ A net (p h ,M h ,Z) M h 

dt XI Mh — M c Ph ’ 


(ID 


where A ne t is the cooling function for gas of metallicity Z. We will 
assume constant temperature, T c — 1000 K, for cold clouds. We 
will use this instantaneous recycling approximation to fix the model 
parameters tg, p t h, esn, and Ao so as to reproduce the Kennicutt 
law in subsection l2.2.3l 


2.2.2 Removing the instantaneous recycling approximation 

We do not use the IRA in our cosmological simulations. In this sec¬ 
tion, we describe how we model feedback and chemical evolution. 
We assume that each stellar particle represents a single stellar pop¬ 
ulation, and so star formation must be treated statistically. Using 
equation j6|, a qualifying particle spawns a new stellar particle of 
mass m * = m 0 iig/N g during a timestep At with probability: 
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We use N g = 3 and have confirmed that our results do not change 
for 7V g > 2. When an SPH particle spawns a stellar particle, the 
mass of cold clouds is reduced by m» and the following new spe¬ 
cific energy and new density p' (in the code we compute the 
effective entropy instead of u e ft) are given to the particle assuming 
the particle volume Msph/p remains constant during the timestep: 


W e ff 


MsPHUeff — m*u c 
Msph — m* 


p = p 


Msph — m* 
Msph 


We then update Msph. 

Each stellar particle has its own age, metallicity, and IMF. We 
define the IMF as the number of stars per logarithmic interval of 
stellar mass per unit total mass of stars formed: 


4>{m) = diV/d In m oc m x \mi<m<m u , 


(13) 


where mi and m u are the lower and upper mass limits of the IMF. 
The IMF has the following normalisation: 



m(f>(m)d\nm 


= 1. 


(14) 


The recycled mass fraction from an evolving population of 
stars is a function of time. For a population of stars formed with an 
IMF 4>(m) at time t = 0, the mass fraction returned to the ISM by 
time t is 


t) 



M r (m)\<j>(m) 


dm 
m ’ 


(15) 


where M(t) is the initial mass of a star just reaching the end 
of its lifetime at t, and M r (m) is the mass of the remnant 
left by _a_star wit h initial mass m . The lifetimes are taken from 
[Portinari. Chiosi, & Bressar] i 19981) for massive stars and iMarisid 
( 2001 ) for intermediate and low mass stars. 

The cumulative number of SNe II explosions up to time t (per 
solar mass of stars formed), Rn{ ^ t), and the mass fraction of 
metals of the j-th element, Ei (^ t), are given by 


t ) 

t) 


rm u 

J max(M(t) ,8Mq ) 





dm 


(16) 

(17) 


where Mi(m) is the mass of the z-th element released from a 
star with mass m. As Ei includes metals present when the stars 
formed, the pre-existing metals must be subtracted in order to esti¬ 
mate chemical yields, Pi(t). If the metallicity of stars with respect 
to the z.-th element at their birth is Z°, Ei can be divided into two 


terms, 



Ei (^ t ) = 

Pi t) + Zi E(^. t), 

(18) 

where 



t) = 

n ( w i dm 

/ yi(m)4>(m) -, 

JM(t) m 

(19) 


and yi (m) indicates the mass of newly produced z-th element in a 
star with mass m. Although Mi and y, depend on the initial stellar 
metallicity in principle, we neglect this dependence and use the 
values for solar initial metallicity throughout, which are taken from 
Portinari et al. (1998). 

We consider energy and chemical feedback from SNe la as 
well. We calculate the number of SNe la explosions using the 


scheme of lGregeio & Renzirfi ti 1983), with parameters updated ac¬ 
cording to Portinari et al. (1998). The progenitors of SNe la are 
assumed to be binary systems with initial masses in the range 
rriB.iow < m B < niB,n P , where ran = mi + m 2 , and mi and 
m 2 are the initial masses of the primary and secondary, respec¬ 
tively. We assume mB, U p = 2mi, ma x, where mi, max is the largest 
single-star mass for which the endpoint is a C-0 white dwarf. The 
binary star systems that are the progenitors of SNe la are assumed 
to have an initial mass function S^>(m B ), where </>(m) is the same 
as for the single-star IMF. The distribution of mass fraction for the 
secondary, p = m 2 /m B , is assumed to have the form (normalised 
over the range 0 < p < 1/2) 

f(p) = 2 1+7 (1 + 7)/z 7 . (20) 


For a single generation of stars formed at t = 0, the number of SNe 
la explosions up to time t is then given by 


/ "IB, up 

</>(m B ) 

l B,low 

where the lower limit 


fl/2 


f(p) ^ 


Mmin (t) 


dm B 

m B 


/Zmin(f) — max 


M(t) m B - m B ,up/2 


m B 


m B 


, ( 21 ) 


( 22 ) 


is set by the conditions that the secondary has evolved off the main 
sequence and that mi = m B — m 2 ^ mi,max. Following Portinari 
et al. (1998), we adopt m B ,iow = 3M 0 , m B , U p = 12Mg, 7 = 2, 
and B = 0.07. Finally, the yield of metals from SNe la is computed 
as 


Pia,i(^ t) = M, Ia _Ri a « f), (23) 

where is taken from the W7 model of lNomoto et alJJl997l) . In 
our code, we follow the evolution of oxygen, iron, magnesium, and 
silicon. Since we will vary the IMF, the self-consistent treatment 
of metals is very important, and allows us to make predictions, for 
example, for the abundance and abundance ratios of the a-elements 
to the iron-peak elements. These aspects of the simulation will be 
considered in a later paper. 

As described above, each stellar particle returns mass, energy, 
and metals to the ambient gas according to its age and metallic¬ 
ity. We smoothly distribute these quantities amongst the neigh¬ 
bouring gas particles using the SPH smoothing. There are vari¬ 
ous possible ways to define the smoothing length of a stellar par¬ 
ticle. As the masses of SPH particles are variable, we choose to 
set the smoothing length of a stellar particle so as to have a con¬ 
stant mass (rather than number) within the smoothing volume. This 
method ensures numerical convergence if one chooses the constant 
mass as M ng b = am OI i g , where M ng b is the gas mass within the 
smoothing volume and a is a parameter whose value should not be 
changed by the resolution. In general, smaller a results in stronger 
feedback effects as expected. As a rational choice we use the same 
mass resolution as in the hydrodynamic calculation, that is a = 40. 
The mass ejected from evolved stars and SNe la is simply added to 
the mass of the hot phases of the neighbouring gas particles. For 
an SPH particle that receives an ejected mass AMfb during the 
timestep, we assign a new effective specific energy u eS and a new 
density p' according to 

1 Msph 

Weft - «eft MgpH - AMpB 

and 

, Msph + AMfb 

p =p 


Msph 
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Using equation the mass evaporated by SNe is calculated 
according to the energy, AQfb. received during the timestep At : 

AM ev = A^TR, (24) 

MSN 

where we compute usn considering only contributions from stars 
heavier than 8Mg using the IMF for the quiescent star formation 
mode (see H2.2.3> . 

The new mass of the hot phase, Mf is given by solving the 
thermal energy equation implicitly, using eqn. (11): 


Mh = Mh + AMev — 


Anet(Phi Mh; z) ^ 

Mh - Me p'h 


(25) 


where the new density of the hot phase, p' h = ph(M^), is given 
by equations f3), J4), and 0 assuming u h and the particle volume 
remain constant during cloud growth. 

The new specific energy of the hot phase, is also calculated 
implicitly, 


M h = Mh + 


AQfb — (tic — mJJAMev 

ML ' 


(26) 


We impose an additional timestep criterion, At fb , for stellar 
particles so as to capture the non-instantaneous feature properly: 


A£ F b = max(cFBf ag e, 0.5Myr), 


(27) 


where cfb is a parameter and Ifb is the age of a stellar population. 
We have confirmed that the simulation results converge for cfb < 
0.1, and thus employ cfb =0.05 throughout. 


Table 1 . The values of the model parameters, adjusted to reproduce the 
Kenni cutt (1998) law represented by eqn (25) in |S 2 nngel_&HernciuisJ 
1200311 . 


Pth i2 

^0 

2 x 10 -25 g cm -3 2.8 Gyr 

1300 



2.2.3 Quiescent star formation mode and parameter setting 

The model described above has three parameters: the threshold 
density, p t h, the evaporation efficiency parameter Ao, and the char¬ 
acteristic star formation timescale t®. Before we can specify these 
parameters, we have to fix the IMF and the amount of energy re¬ 
leased by a SN. For this quiescent mode of_star formation, we 
use the Salpeter IMF (x = 1.35 in eq. < 1 3t : ISalnetertl 195.4l with 
m u = 100 Mq and mi = O.IMq and adopt 2 x 10 51 erg for 
the energy released per SN, which yields a SN temperature Tsn = 
2pustt / (3k) ~ 3.3 x 10 8 K, where p denotes the mean molecular 
weight. The reason why we use this rather than the canonical value 
of 10 S1 erg per SN is simply that we need this amount of energy 
in order to reproduce the Kennicutt (1998) law for gas with solar 
metallicity. For primordial gas, the value 10 51 erg/SN is sufficient. 
Note that an IMF with a single slope is too simple a description 
Ie.g. lKennicuttll9 83: Kroupa 1998) and changing the slope of the 
low mass portion, say for in < 1 Mg, can easily change the number 
of SNe II by a factor of ~ 2. For simplicity, we restrict ourselves 
to using a single-slope IMF, but we reserve the freedom to increase 
the energy per supemovae just discussed. 

To fix the model parameters, we considered self-regulated star 
formation in a self-gravitating gas sheet, adopting the instantaneous 
recycling approximati on. This procedure is identical to that fol¬ 
lowed by Springel & Hernauist (2003), so we will give no further 
details here. The values of the parameters we adopt are given in 
Table Q 

In Fig.0 we show the relation between the star formation sur¬ 
face density and the gas surface density in an idealis ed sim ulation 
of disc formation in a virialised halo (see lOkamoto et all 1200.1 ). 
Even though the values of the parameters were set assuming the 
instantaneous recycling approximation and solar metallicity, this 
simulation, which did not make this approximation and followed 


Figure 1. Star formation rate per unit area versus gas surface density in 
an idealised simulation of disc formation in a virialised halo. The surface 
star formation rate densities are computed using the surface density of stars 
younger than 3 X 10 7 yr in cylindrical bins. The outermost bin corresponds 
to the edge of the star forming region. The symbols represent star formation 
rates at t = 0.25, 1.5, 2.75, 4, 5.25 Gyr, indicated by crosses, diamonds, 
triangles, squares, and crosses respectively. The solid and dotted lines indi¬ 
cate the target relation and the cut-off surface density respectively. 

chemical evolution, shows reasonable agreement with the target re¬ 
lation. The surface gas density at the edge of star-forming region is 
also consistent with the threshold surface density. 

2.2.4 The burst mode and top-heavy IMF 

In the self-regulated star formation that occurs in the quiescent star 
formation mode, cooling and feedback are balanced. We implement 
a second mode of star formation, which is able to blow gas out 
of galaxies. For this to happen, the injection of feedback energy 
has to be more rapid. Hence, we require that star formation should 
occur in short bursts in physical conditions that permit heating by 
feedback to exceed local cooling. After various tests, we found that 
sufficient heating cannot be obtained simply by adopting a shorter 
t* 0 (or equivalently, a larger c„). Motivated by the considerations 
discussed in § 1, we therefore decided to assume a top-heavy IMF 
in the burst mode as well as a short star formation timescale. We 
use to — 0.15 Gyr which corresponds to c* = 0.5 and an IMF 
with slope x = 0.34, which maximises the number of SNe II. The 
other parameters including the upper and lower mass limits for the 
IMF are the same as for quiescent star formation. 

In semi-analytic galaxy formation models, it is usually as- 
sum ed that a starburst is triggered by a major galaxy m erger 

ie.g. lWhite & Frenkll99ll : !Kauffmann. White & Guiderdonij 199.31: 
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Cole et a il boocb , although some studies include small starbursts 
induced by minor mergers as well ISomerville & Prim acid 1 1 99^ : 
lOkamoto & Napashirnall2003ll . In our simulations, we would like 
to model the burst on the basis of local physical quantities rather 
than on global information such as the merger mass ratio, which is 
not easily available. Exactly how to do this, however, is not evident. 
In this paper, we examine two possibilities for triggering the burst 
mode of star formation. 

_ Firstly, we consider the gas density. iMihos & Hernauistl 

1 19941) showed that cold gas is driven to the centre of the rem¬ 
nant when galaxies merge. The burst is therefore preceded by an 
increase in the central density and it seems plausible that imposing 
a threshold density, pburst, above which the burst mode is switched 
on, can capture the conditions under whic h nuclear starbursts oc- 
cur in merging galaxies. In addition. ISpringel & Hernauistl 120031) 
noted that self-regulation in the multi-phase model breaks down at 
sufficiently high density; our “density-induced” burst may be con¬ 
sidered as a crude modelling of this process. 

The second plausible trigger for bursts is the presence of 
shocks. Galaxy mergers induce strong, galactic-scale shocks. In¬ 
deed, iBarnesI l2004t) suggested that shock-induced star formation 
is the dominant mode in interacting galaxies. We model a shock- 
induced burst as follows. When the rate of change of the entropy 
variable K(s ) = (7 — 1 )up 1-7 due to the artificial viscosity ex¬ 
ceeds a threshold value, I\ burst, the burst mode is switched on. 
We use 7 = 5/3 throughout. Note that other variables such as u 
or p/(p fdyn) can also be used to identify shocked particles, and 
we have confirmed that these alternative choices produce the same 
results for appropriate values of the thresholds. In contrast to the 
density-induced burst, the shock-induced burst is expected to result 
in a large scale (galactic scale) starburst in interacting galaxies. The 
shock trigger is more sensitive to mergers than the density trigger 
because in the latter case, the burst cannot start until sufficient gas 
has been funneled into the galactic centre by the merger dBarnesI 
•:2004). 

We use pburst = 10 -23 g cm” 3 and TTburst = 9 x 10 3 in 
our code units ( h _1 Mpc, km s _1 , and 10 1 o /i _1 Mq for length, 
velocity, and mass, respectively) throughout this paper. The depen¬ 
dence of our results on these parameters will be discussed later. It 
should be noted that, for the density-induced burst, our results do 
not depend on resolution as long as this is good enough to resolve 
the threshold density pb urs t- On the other hand, the shock-induced 
burst will have a strong dependence on resolution because the width 
of the shocked layer is simply proportional to the smoothing length. 
Thus, we would need to adjust Tvburst or tg for the burst mode if 
we were to change the numerical resolution. 


3 SIMULATION SETUP 

In order to study disc formation, we have selected from a pre¬ 
existing cosmological IV-body simulation a halo with_a quiet 
merger history. The semi-analytic model of lCole et alJ ( 200d ) ap¬ 
plied to the merger tree of this halo, predicts a galaxy that is disc- 
dominated at the present d ay. This halo is, in fact, the same as the 
one used in lOkamoto et alJ 1200 3). In that paper, we showed that, if 
cooling is not allowed to occur until 2 = 1 , then a reasonable disc 
galaxy forms by 2 = 0. 


3.1 Initial condition 

We assume a low-density, flat CDM universe (ACDM) as the back¬ 
ground cosmology. We adopt the following choices for the cosmo¬ 
logical parameters: mean matter density, Ho = 0.3, Hubble param¬ 
eter, h = 77o/100kms _1 Mpc _1 = 0.7, cosmological constant 
term, Ha = Ao/(3 Hq) = 0.7, amplitude of mass fluctuations, 
<78 = 0.9, and mean baryon density, Hb = 0.04 

To generate our initial^conclitions, we use the resimulation 
technique introduced by iFrenk et all (.1996). The halo of interest 
was grown in a dark matter simulation of a 35.325 /i” 1 Mpc peri¬ 
odic cube. Its mass is about M v n = 1.2 x 10 12 /i _1 Mq within the 
sphere which has the virial overdensity, 5 v i r = 337, at z = 0. The 
circular velocity, spin parameter, and collapse redshift of this halo 
are v c (r vil ) = 155 km s” 1 , A = J|£j 1/2 /(GM 5/2 ) = 0.038 
and 2 C ~ 1.5, respectively, where z c is defined as the redshift at 
which the main progenitor had half the final halo mass. To gener¬ 
ate the new initial conditions, the initial density field of the par¬ 
ent simulation is recreated and appropriate additional short wave¬ 
length perturbations are added in the region out of which the halo 
forms. In this region we also place SPH particles in the ratio of 1:1 
with dark matter particles. The region external to this was popu¬ 
lated with high-mass dark matter particles, the function of which is 
to reproduce the appropriate tidal fields. The initial redshift of the 
simulation is 50. The masses of the SPH and high-resolution dark 
matter particles are ~ 2.6 x 10 6 h _1 Mq and ~ 1.7 x 10 7 /i” 1 Mq, 
respectively. 

The gravitational softening lengths are kept fixed in comoving 
coordinates for 2 > 3; thereafter they are frozen in physical units 
at a value (of the equivalent Plummer softening) of t = 0.5 kpc 
and 1 kpc for the SPH and high-resolution dark matter particles 
respectively. The gravitational force obeys the exact r” 2 law at r > 
2.8e. 

3.2 Models 

To investigate the effects of the burst mode of star formation on 
the evolution of the galaxy, we consider three models. In the first, 
we do not include a burst mode. We refer to this reference case as 
the no-burst model. The second (the density burst model) has the 
density-induced burst mode and the third (the shock-burst model), 
the shock-induced burst mode. To study the effects of the burst 
mode in detail, we also consider, in 44.41 models with different 
threshold values for triggering the burst, as well as models with 
a combination of density-induced and shock-induced bursts. 


4 RESULTS 

First, we present a description of the final galaxy that forms in each 
simulation with a different feedback prescription, and then under¬ 
take a detailed study of the evolution in each case. In the final sub¬ 
section, the burst parameters are varied in order to test the sensitiv¬ 
ity of the simulation to these assumptions. 

4.1 Galaxy properties at 2 = 0 

Fig- HI shows the stellar and gas distributions within 50 /i” 1 kpc 
boxes centred on the galaxies at 2 = 0. The edge-on and face- 
on projections are selected to be perpendicular and parallel to the 
angular momentum vector of the stellar component within the cen¬ 
tral 10 hr 1 kpc sphere. The no-burst galaxy has a tiny stellar disc 
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Figure 2. Galaxies at z = 0. The no-burst, density-burst, and shock-burst models are shown from left to right. The edge-on and face-on views of stars, and 
edge-on and face-on views of gas are given in rows from top to bottom. We use the nett angular momentum of stars within 10 h _1 kpc spheres to define the 
viewing angles. The brightness indicates the projected mass density, and the same scaling is used for each model. All of these galaxies are obtained from the 
same initial conditions. 
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Figure 3. Upper panel: Surface brightness profiles in the /-band for the no¬ 
burst (crosses), density-burst (diamonds), and shock-burst (triangles) mod¬ 
els. The solid lines are double-exponential fits for r < 20 kpc and the fit¬ 
ting parameters are given in Tabic [/l Lower panel: Cold gas surface density 
profiles. The solid, dotted, and dashed lines represent the no-burst, density- 
burst, and shock-burst models, respectively. Note that we have removed the 
h dependence here. 

and a large population of halo stars. Visually, this galaxy is similar 
to the one that lAbadi et alJ (;2003a) obtained in their simulation. It 
should be noted that our effective equation of state for the ISM is as 
stiff as that used by R obertson et alJ 120041) . and hence the stabili¬ 
sation of gas against the Toomre instability is insufficient to allow 
a large stellar disc to form in this particular halo. The simulation 
with the density-burst model produces an ellipsoidal stellar object 
that is slowly rotating. Due to strong feedback in the galactic centre 
arising from the density-induced starbursts, most of the gas is ex¬ 
pelled from the galaxy and only a diffuse gaseous ring remains. In 
contrast, in the shock-burst model the number of halo stars is sig¬ 
nificantly lower and an extended disc forms in both the stellar and 
gaseous components. At the galactic centre, dense gas cores always 
exist regardless of the model. We suspect that this is an artificial ef¬ 
fect caused by our limited resolution. 

The /-band surface brightness profile of each galaxy is shown 
in the upper panel of Fig. [3] To compute luminositieSjjve used the 
popul ation synthesis code PEGASE2 jFioc & Rocca-Volmerangel 
119971) . The luminosity of each particle is calculated for the IMF, 
metallicity and age appropriate to that particle in the simulation. 
We fit each surface brightness profile up to 20 kpc with a double- 


Table 2. Values of the parameters used in Fig. to fit the /-band sur¬ 
face brightness profiles. The central surface brightnesses S° and the scale 
lengths R d are given for the inner and outer profiles. Subscripts i and o in¬ 
dicate the inner and outer profiles respectively. The last 2 columns give the 
total B- and /-band luminosities for each galaxy. 



s? 

M 

0 o 

Rf/kpc 

-R^/kpc 

M b 

Mi 

no-burst 

19.1 

18.8 

1.5 

3.3 

-21.1 

-22.6 

density-burst 

19.1 

20.6 

0.8 

4.0 

-19.7 

-21.4 

shock-burst 

17.7 

22.1 

1.0 

39 

-21.3 

-22.7 


exponential because the standard r 1 ^ 4 + exponential profile cannot 
fit the shock-burst galaxy. It should be noted, however, that the r 1//4 
+ exponential profile is a better description for the density-burst 
model, and therefore the double-exponential fit underestimates its 
bulge. Values of the fitting parameters are shown in Table[2]which 
also gives the total B- and /-band luminosities of each galaxy. 

The surface brightness profile of the no-burst galaxy in Fig. [3] 
is well fit by a single exponential (90 % of the /-band light comes 
from the outer exponential disc). Thus, there is no evidence in this 
profile for a bulge component even though Fig. [2] clearly shows 
an extended spheroidal component around this galaxy. The double¬ 
exponential is a good fit to the shock-burst galaxy, but the outer 
exponential is rather flat, indicating that the surface brightness of 
the disc is nearly constant as a function of radius. The resulting 
scale length, 39 kpc, is therefore very large. As we discuss in §5, 
it is possible that the disc acquires excessive angular momentum 
because the galactic winds are not collimated and exert a pressure 
on the hot gas reservoir which becomes distended and more sus¬ 
ceptible to tidal torques at early times. In spite of its large size, the 
shock-burst galaxy falls on the /-band Tully-Fisher relation. For its 
rotation velocity of v TO t — 200 km s -1 , its absolute magnitude, 
Mi = —22.7, is only slightly brighter than the mean relation and 
well within the scatter iGiovanelli et al . 1991. 

The surface brightness profile in each galaxy reflects the cold 
gas surface density profile (lower panel). In the no-burst galaxy, 
the cold gas disc has an approximately exponential profile that is 
slightly more extended than the distribution of the stellar light. The 
density-burst galaxy has lost its gas content due to strong winds 
as we will show in §4.2. As a result, the gas surface density is an 
order of magnitude lower than the Kennicutt threshold (T, gas ~ 
10 Mq pc -2 ). In the shock-burst model, a plateau is also seen in 
the gas surface density profile, and the surface density is close to 
the Kennicutt threshold. The three different models all have similar 
central surface gas densities. This implies that the central gas cores 
are a numerical artifact and higher resolution or kinetic feedback 
may be required to get rid of these cores. 

The fits to the surface brightness profiles of Fig.|3|provide only 
partial information about the morphology of the galaxy. A more in¬ 
formative way to characterise the relative importance of bulge and 
disc components is to carry o ut the dynamical decomposition pro¬ 
posed bv lAbadi et alJ<2003bl) . For this, we first compute the angu¬ 
lar momentum, J z , of each star particle parallel to the nett angular 
momentum of stars within 10/i -1 kpc, and the angular momentum 
of the co-rotating circular orbit, J C (E). The ratio J z /J c defines 
an orbital circularity. In Fig0 we show the probability distribu¬ 
tion of this orbital circularity for stars within 25 fi -1 kpc of the 
galactic centre. A disc component should have J Z /J C (E) ~ 1 and 
such a component is clearly visible in the shock-burst model, con¬ 
firming the visual impression gleaned from Fig. [2] By assuming a 
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Figure 4. Mass-weighted probability distributions of the orbital circularity, 
J Z /ME), within 25 h ~ 1 kpc from the galactic centre. The solid, dotted, 
and dashed lines indicate the no-burst, density-burst, and shock-burst mod¬ 
els respectively. 


Table 3. Disc-to-total mass and luminosity ratios for the simulated galaxies 
at z = 0. 


5 3 2 1.5 1 Z 0.5 0 




mass 

U 

B 

V 

I 

K 

no-burst 

0.26 

0.72 

0.63 

0.54 

0.45 

0.44 

density-burst 

0.21 

0.22 

0.22 

0.22 

0.21 

0.20 

shock-burst 

0.48 

0.86 

0.84 

0.80 

0.72 

0.66 


Figure 5. Upper panel: the formation history (star formation rate as func¬ 
tion of time) of the stars that lie within 25 h~ 1 kpc from the galactic centre 
at z = 0. Lower panel: the fraction of stars formed in the burst mode. The 
solid, dotted, and dashed lines correspond to the non-burst, density-burst, 
and shock-burst models respectively. 


non-rotating spheroid, i.e. that stars in the spheroid are symmet¬ 
rically distributed around zero, all stars having J Z /J C {E) Sj 0 
are identified as a half of the spheroid. Their counterparts with 
Jz/Jc(E) > 0 are defined statistically. All remaining stars are 
identified as the disc component. We do not try to decompose the 
disc into thick and thin discs but defer a detailed study to a forth¬ 
coming paper. 

In Table. 0 we show the disc-to-total ratios for mass, U, B, 
V, I, and A'-bands. The redder the band, the more the spheroid 
dominates. These ratios confirm that the shock-burst model has the 
most significant disc, while the density-burst model is the most 
spheroid-dominated. The 5-band disc-to-total luminosity ratio of 
the shock-burst galaxy is D/T = 0.84, which is sufficiently large 
to be identified conventionally as a disc galaxy, while the rela¬ 
tion betwee n the D/T s and the Hubble T-ty pes has a large scatter 
(Baugh. Cole. & FrenlJl99fiHGrahaml200ll) . Note that the dynami¬ 
cal decomposition reveals a significant spheroidal component in the 
no-burst galaxy even though the photometric decomposition based 
on the surface brightness profile failed to detect it. 

4.2 Mass accretion and star-formation histories 

The star formation histories (SFHs) of the simulated galaxies are 
plotted in Fig.|5] In the no-burst simulation, significant star forma¬ 
tion activity takes place at high-redshift, and this early star forma¬ 
tion is responsible for the massive bulge (and halo) population that 


forms in this model. This galaxy has a nearly constant star forma¬ 
tion rate after z — 1, and it builds up the tiny disc present at the 
final time. 

The SFH of the density-burst model is similar to that of the 
no-burst model up to the point when the gas density reaches the 
threshold for a burst. At this point, the burst is triggered, the asso¬ 
ciated feedback suppresses further star formation, and the star for¬ 
mation rate falls well below that in the no-burst model. The huge 
amount of feedback energy released into a small region near the 
centre of the galaxy blows out most of the hot gas and this al¬ 
most completely quenches further star formation. As a result, the 
density-induced burst simulation yields an elliptical-like galaxy. If 
the threshold is lowered, then the burst is turned on earlier and even 
fewer stars form. Thus, it is difficult to produce a galaxy that is 
more disc-dominated than in the no-burst case simply by allowing 
density-induced bursts. 

By contrast, shock-induced bursts occur even at very high red- 
shift, in small haloes where bursts are triggered both by mergers 
and violent collapse. The strong feedback from the stars formed in 
these bursts suppresses cooling and star formation during this early 
period when the galaxy is undergoing a number of merger events. 
As a result, the peak in the SFH shifts to t ~ 5.3 Gyr (z ~ 1.1). 
An extended reservoir of hot gas is created by the feedback energy 
associated with the burst and much of it remains attached to the 
halo. After z = 1, when the bulk of the dynamical activity is over, 
this gas cools to form an extended , young stellar disc. The key fea¬ 
tures of the shock-induced burst may be seen in the lower panel of 
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Figure 6. Evolution of the mass in the main progenitor haloes. The solid, 
dotted, and dashed lines represent the normalised halo mass, baryon mass, 
and the cold baryon (cold gas and stars) mass in the halo, respectively. We 
show the no-burst, density-burst, and shock-burst models from top to bot¬ 
tom. 

Fig- El The fraction of stars born in the burst mode is large (~ 40 
per cent) at high-redshifts (t < 6 Gyr) and it gradually decreases 
to ~ 5 per cent towards the present. Thus, the shock-burst model 
experiences strong feedback at early times, followed more recently 
by a sustained period of steady, self-regulated star formation. 

The top-heavy IMF present in bursts delivers feedback en¬ 
ergy in a particularly efficient way. It is worth noting that even at 
its peak, the star formation rate in the shock-burst model barely 
reaches ~ 20 Mq /yr. However, the star formation rate in the star- 
burst would be significantly overestimated if, as is the norm in ob¬ 
servational studies, the rates are inferred by assuming a standard 
IMF. 

To investigate the build-up of the visible components of the 
galaxy and the importance of galactic winds, we now consider the 
mass evolution of the main progenitor of the halo as a function 
of time. We define the mass of the main progenitor as Mtot = 
(47r/3)r®j r /9 V ir(z) where r v i r i s given by the spherical c ollapse 
model in our chosen cosmology lEke. Cole & Frenkl ^96). Since 
the accretion and cooling of baryonic matter affects the structure 
of the dark matter halo, the halo mass defined in this way differs 
slightly between models: the more concentrated the dark matter is, 
the smaller the halo mass. In Fig. [6] we plot the normalised main 
progenitor halo mass, /bMh a io/(l — /b), where /b = fib/fio is 
the baryon fraction in the universe. Also plotted are the mass of 
baryonic matter, Mb, and the mass of cold baryons (stars and cold 
gas), Mcoia, in the main progenitor halo. If the halo had the mean 
cosmic baryon fraction, /b, then the normalised halo mass would 
be equal to the baryon mass. 


Table 4. The stellar mass formed within a comoving 25// —1 kpc sphere 
around the centre of the main progenitor divided by the total stellar mass 
within 25/i _1 kpc from the galactic centre at z = 0. 


no-burst 

density-burst 

shock-burst 

0.60 

0.25 

0.78 


The baryon fraction in the no-burst galaxy is always very close 
to the cosmic value. Very little gas is lost from the halo in this 
simulation because all the feedback energy goes mostly to sup¬ 
port the ISM against its self-gravity. In the density-burst simulation, 
the baryon fraction in the halo is quite low and there is almost no 
hot gas left, signaling the existence of strong galactic winds. The 
stepwise jumps in the cold baryon mass indicates that this mass 
grows primarily through mergers with other haloes. The shock- 
burst galaxy also has winds, but these are not as strong as in the 
density-burst model. By the final time, the halo of the shock-burst 
galaxy has lost 24 per cent of its baryons. While the cold baryon 
content at z = 0 is similar to that in the no-burst simulation, the 
rate of growth of the cold baryons is steeper at early times in the 
no-burst model and at late times in the shock-burst model. This is 
the origin of the difference in the morphology of the final galaxy in 
these two cases. 

The behaviour of both the SFH and the mass evolution in the 
shock-burst model indicates that the success of this model stems 
from the suppression of star formation in small haloes. In order to 
identify the location of the main star formation sites, we plot, in 
Fig. El the birthplaces, relative to the position of the main progen¬ 
itor at the time, of the stars that end up within 25 hr 1 kpc of the 
galactic centre at z = 0. The continuous trajectories of the birth¬ 
places in the no-burst model (top left panel) indicate that small 
haloes undergo continuous star formation until they fall into the 
main progenitor. A similar behaviour is seen in the density-burst 
model. Since density-induced bursts do not occur in small haloes, 
where the gas density cannot reach the threshold value, star forma¬ 
tion in such small systems is the same as in the no-burst model. The 
main difference between the two models is apparent in the zoomed 
panel on the right-hand side of the figure: strong feedback from 
density-induced bursts almost stops star formation in the main pro¬ 
genitor after 2 ~ 1. The main sites of star formation in this model 
are small infalling haloes. This is how the elliptical galaxy forms 
in the density-burst model. Unlike density-induced bursts, shock- 
induced bursts can occur in small haloes whose shallow potential 
wells make feedback particularly effective. The star formation tra¬ 
jectories are therefore no longer continuous. Instead, intermittent 
star formation in infalling haloes is clearly evident in the lower 
left-hand panel. As anticipated, the suppression of star formation 
in small galactic building blocks is the key to making a large disc. 

More quantitative information on the location of the star for¬ 
mation sites is provided in Table El which lists the fraction of the 
final stellar mass within 25h~ 1 kpc from the galactic centre that 
has formed within a comoving sphere of radius 25 h^ 1 kpc around 
the main progenitor. The mass of each star used in this calculation 
is the mass at z = 0, not the initial mass. We find that 78 per cent 
of the stellar mass formed in the main progenitor in the shock-burst 
model. In contrast, this fraction drops to 25 per cent in the case of 
the density-burst model. 
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Figure 7. Birthplaces of stars that lie within 25/i -1 kpc of the z = 0 galactic centre. These birthplaces are plotted relative to the position of the centre of the 
main progenitor at the epoch when each star forms. The colours indicate the formation time of the stars, with red being the oldest and white the youngest. Panels 
show the no-burst (top row), density-burst (middle), and shock-burst (bottom) models. The left and right panels show the x-y projection of the birthplaces in 
2h~ 1 Mpc and 200/i —1 kpc cubes respectively. The coordinate system used here is that of the original simulation, i.e. it has not been rotated to match the 
angular momentum of the final galaxy. 
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Figure 8. Gas distribution in physical 50 h 1 kpc boxes centred on the main progenitors in the shock-burst model. A fixed viewing angle is used, that for 
which the disc is edge-on at z = 0. The z = 0 gas distribution is plotted in the bottom-right panel of Fig. [3 
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Figure 9. Evolution of the angles between the angular momentum vec¬ 
tor of the stellar component within 10h~ 1 kpc at z = 0 and the angular 
momentum vectors of various components in the main progenitor for the 
shock-burst simulation. The different lines refer to the following compo¬ 
nents of the main progenitors: stars within 25/r~ 1 kpc (solid), cold gas 
within 25 h ~ 1 kpc (dotted), dark matter within 25 h ~ 1 kpc (dashed), and 
dark matter within the virial radius (dot-dashed). Angles are divided by 7r, 
and 9 /tr = 1 means counter-rotation relative to the final disc. 

4.3 The evolution of angular momentum 

Since the direction of the angular momentum vector of the main 
progenitor halo varies in time due to gravitational torques, it is ex¬ 
pected that the direction of the angular momentum vector of the 
accreting gas will also vary. Such a process may affect the forma¬ 
tion of the disc through angular momentum mixing. In Fig. [8] we 
show the gas distribution around the main progenitor in the shock- 
burst model. The viewing angle is fixed and chosen so that the disc 
is edge-on at 2 = 0 (see the bottom-right panel of Fig. [3. A gas 
disc already exists at 2 = 1, but its orientation changes signifi¬ 
cantly throughout its evolution. For example, the disc is inclined 
and almost edge-on at 2 = 1, but it becomes almost face-on at 
2 = 0.7 before settling down to the final orientation at 2 = 0.2. 
A good example of angular momentum mixing is seen in the panel 
for 2 = 0.4, where the newly accreting gas does not line up with 
the pre-existing inner disc but settles instead onto a different plane 
and the two discs torque one another. As these gas discs host star 
formation, this process contributes to the formation of hotter stellar 
components, such as a bulge or a thick disc. 

In Fig. [9] we plot the evolution of the direction of the angu¬ 
lar momentum vectors of various components of the main progen¬ 
itor in the shock-burst model. The lines show the angles between 
each angular momentum vector and that of the final stellar disc. 
Angular momenta are calculated for material inside a sphere of a 
radius 25/i _1 kpc (in physical coordinates) for stars (solid), cold 
gas (dotted), and dark matter in the core (dashed), while the angu¬ 
lar momentum of the main progenitor halo is defined using the dark 
matter inside the virial radius (dot-dashed). At early times, t < 8.5 
Gyr, the spin of the cold gas correlates well with that of the dark 
halo, while the stars show a better correlation with the core at t < 6 
Gyr. Note that the angular momentum of the stellar component is 
contaminated by satellites and hence it is rather noisy. On the other 
hand, at late times, t > 8.5 Gyr, corresponding to 2 < 0.5, the 
stars follow the cold gas quite well and they are offset from both 
the dark matter core and halo. 


We confirmed that the angular momentum flip of the cold gas 
disc seen in Figs. |8]is caused by minor mergers occurring between 
2 ~ 1-0.6. Afterwords, its evolution is driven by angular momen¬ 
tum mixing with newly accreted gas that has a different orientation 
from the pre-existing disc (see Fig. |8j. For the stellar component, 
initially the direction of its angular momentum vector is well cor¬ 
related with that of the dark matter in the core, both of which result 
from merging and dynamical friction. As stars begin to form in 
the well-developed gas disc, the nett stellar angular momentum be¬ 
comes dominated by disc stars, and the angular momentum vectors 
of the stellar and cold gas component become aligned. 

4.4 Dependence on burst parameters 

In this subsection, we discuss how galaxy properties are affected 
by changes in the burst parameters. For the shock-burst model, we 
have checked that if we adopt the same star formation efficiency in 
the burst mode as in the quiescent mode, but retain the top-heavy 
IMF, then the results are very similar to the no-burst simulation. 
Thus, a top-heavy burst IMF alone is insufficient to produce a big¬ 
ger disc. Similarly, merely increasing the star formation efficiency 
in the burst without including a top-heavy IMF does not lead to a 
significantly bigger disc than is found in the no-burst simulation. 
Only the combination of these two features, i.e. a top-heavy burst 
IMF and an increased burst star formation efficiency (c* ~ 1 in 
eqn. 2) is capable of promoting stronger feedback to the extent that 
an extended disc can form in this particular halo. 

Since the response of the density-burst model to changes in 
the density threshold is fairly straightforward and the shock-burst 
model looks more promising, we decided to explore the following 
three additional models: 

• SFI: The same as the shock-burst model, but with a slightly 
higher threshold for the burst, A'burst = 10 4 , rather than 7\ burst = 
9 x 10 3 as in the standard shock-burst model. 

• SL: The same as the shock-burst model, but with a lower 
threshold for the burst, A'burst = 8 x 10 3 . 

• DS: A combination of density-induced and shock-induced 
burst modes, with the same parameters as the original models. 

In Fig. 1101 we show the SFHs of the galaxies in these three 
variant models. Models with shock-induced bursts are quite sen¬ 
sitive to the threshold, Aburst- In general, the lower the thresh¬ 
old, the fewer the number of stars that form because of the ad¬ 
ditional amount of feedback generated by the burst stars. The 
DS model shows interesting behaviour. As the feedback from the 
shock-induced bursts suppresses gas cooling at early times, the red- 
shift at which the gas density reaches the threshold for a density- 
induced burst is reduced. At this time, the depth of the potential 
well is deeper and this decreases the strength of galactic winds. As 
a result, when the density-induced bursts start, the star formation 
rate exceeds that in the shock-burst model. These density-induced 
bursts in the DS model do not lead to a smaller final bulge mass, but 
do suppress further cooling, showing how difficult it is to promote 
disc formation when a large amount of feedback energy is pumped 
into the galactic centre after 2 ^ 2. 

Edge-on views of the final galaxies are shown in Fig. IIII 
Shock-induced bursts assist disc formation and lowering the thresh¬ 
old for such bursts results in a larger disc. Of course, too low a 
threshold is counterproductive. For example, for A'burst = 5 x 10 3 , 
the star formation becomes dominated by the burst mode at all 
times. As a result, only a very dim object forms because the thresh¬ 
old is below the typical values in the accretion shock and most of 
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Figure 11. Edge-on projections of stars (upper row) and gas (lower row) in the SH, SL. and DS models at 2 = 0. 


Table 5. Disc-to-total ratios as in Tabic HI only this time for galaxies in the 
SL, SH, and DS models. 



mass 

U 

B 

V 

I 

K 

SH 

0.39 

0.83 

0.80 

0.76 

0.66 

0.59 

SL 

0.55 

0.94 

0.91 

0.87 

0.79 

0.73 

DS 

0.20 

0.51 

0.47 

0.41 

0.33 

0.28 


the stars are then formed outside the central disc. The galaxy in 
the DS model has an outer gas ring. One might suspect that this is 
just gas swept out by feedback. In fact, the ring is newly accreted 
gas that had a different angular momentum direction relative to that 
of the inner disc. When the outer ring is either spatially separated 
from the inner disc or perpendicular to it, it tends to be stable. From 
our experience, such outer rings are not rare in cosmological simu¬ 
lations because the direction of the angular momentum of the main 
progenitor halo can change dramatically, as we showed in Fig. [9] 
and so the orientation of newly accreting gas can differ from that of 
the pre-existing disc (see Fig. [8}. 

We also performed a dynamical decomposition for these 
galaxies (see Table [5j. As expected, lowering the threshold for a 
shock-induced burst leads to larger disc-to-total ratios. These num¬ 
bers also confirm that density-induced bursts result in early type 
galaxies even when they are combined with shock-induced bursts. 


5 SUMMARY AND DISCUSSION 

We have performed cosmological simulations of galaxy formation 
in a cold dark matter halo chosen to have a quiet merger history. 
We have employed a multi-phase description of the ISM based on 
that of lSpringel & Hemauis t: 42003 ) and developed explicit models 
for star formation in quiescent and burst modes. We assume that 
the star formation efficiency is much higher in the burst mode than 
in the quiescent mode and, motivated by the semi-analytic work 
of iBaugh et al. (2005), that the stars in a burst form with a top- 
heavy IMF. Our model potentially provides a unified picture of star 
formation in galaxies of all morphological types. In this paper, we 
have explored two triggers for the starburst mode: strong shocks 
and a high density. An important conclusion of our work is that, for 
the same initial conditions, different burst models result in galaxies 
with a wide range of morphological types, having bulge-to-total 13- 
band luminosities spanning the range D/T ~ 0.2 to 0.9. 

Although, having simulated galaxy formation in only one 
halo, we must be cautious about generalising our conclusions, the 
shock-induced burst model does look promising for creating disc- 
dominated galaxies in haloes with relatively quiet recent merger 
histories. The lower the threshold for shock-induced bursts, the 
more the galaxy becomes disc-dominated. The key to the success of 
this model in a universe where structure grows hierarchically is that 
the burst fraction is high at early times and subsequently decreases, 
once the main phase of halo merging activity has taken place. Feed¬ 
back from the top-heavy IMF in the shock-induced bursts sup¬ 
presses the early collapse of baryons in small haloes, helping to 
create a reservoir of hot gas that remains attached to the halo and is 
available for cooling later on. It is this gas that eventually ends up 
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Figure 10. The same as Fig. [s] but for the SH (solid), SL (dotted), and 
DS (dashed) models. The standard shock-density model is also shown for 
comparison (dot-dashed). 


forming a large, young stellar disc in which stars continue to form 
at a slow rate up to the present day. 

In spite of its success in generating a disc-dominated galaxy 
which falls on the /-band Tully-Fisher relation, our shock-induced 
burst model produced an excessively large disc with an unrealisti¬ 
cally flat surface density profile. A possible cause of this problem 
may be the lack of collimation of the galactic winds in the simula¬ 
tion. The pressure generated by the roughly isotropic winds causes 
the hot gas reservoir to becomes more extended than the dark mat¬ 
ter distribution at early times and thus susceptible to gaining extra 
angular momentum from tidal torques. When this gas later cools, 
it settles onto a disc with too much angular momentum. Thus, far 
from having too little angular momentum as most previously sim¬ 
ulated discs, our disc ended up with too much angular momentum! 
If this interpretation is correct, the lack of collimation in the winds 
might just reflect the inability of SPH to represent the steep pres¬ 
sure gradients required to collimate the outflow. We intend to test 
this hypothesis with simulations using a mesh-based hydrodynamic 
code. This result cautions about the general assumption that the 
baryonic ma tter has the same specifi c angular momentum of the 
dark matter jFall & EfstathioiiFl980 ). In fact, feedback is able to 
change the distribution of the gas in proto-galactic regions, and 
therefore baryons can have higher angular momentum than dark 
matter. 

While other improvements implemented in our simulations, 
such as the phase-decoupling method to suppress numerical trans¬ 


fer of angular momentum, high resolution, and the use of a stiff 
equation of state f or the ISM, also contribute to the form ation 
of large discs (se e lOka moto et alJ 120031 : iGovernato et alJ 120041 : 
iRobertson et al Jl/OOdl) . the comparison between models with and 
without shock-induced bursts shows that our burst model is the key 
ingredient that helps resolve the angular momentum problem in a 
CDM universe. Although the assumption of a top-heavy IMF in 
bursts is controversial, it is encouraging that the model we have 
used is the same that was originally introduced in order to account 
for the abundance an d pro perties of galaxy populations at high- 
redshift feaugh et d Eooa . and was subsequently found also to 
provide a plausible explanatio n for the metal content in the intra- 
cluster medium Nagashimaet air2005^ and in elliptical galaxies 
iNagashima et al]l2005bl)~ 

By contrast, the density-induced burst model seems to make 
things worse for disc formation. Since the gas density only reaches 
the threshold for bursts in massive haloes, this model does not af¬ 
fect gas cooling and star formation in small haloes. The feedback 
associated with the burst merely halts star formation in the galaxy 
at late times when an extended stellar disc might otherwise be ex¬ 
pected to form. Consequently, the density-induced burst leads to 
the formation of an early-type galaxy. The effects of bursts (with 
their top-heavy IMF) in this model may resemble those generated 
by energy fed back from active galactic nuclei (AGN). While AGN 
do not produce stars and metals, both forms of feedback inject large 
amounts of energy into the galactic centre. Thus, we suspect that it 
is difficult to solve the angular momentum problem in simulations 
of disc galaxies solely by including AGN feedback. Of course, it is 
entirely possible that AGN feedback beha ves q uite differently from 
the case of density-induced bursts (se e iKawata & Gibson ti2004) 
and lDi Matteo. Springel & Hernauistl l2005l) for attempts at mod¬ 
elling AGN feedback on galactic scales). However, if feedback 
from AGN does have similar consequences to those produced by 
the density bursts in our simulations, then one might expect that 
it is particularly important for suppressing recent star formation in 
the most massive galaxies. This is in qualitative agreement with the 
observation that the most massive galaxies tend to be ellipsoidal 
and not discy. 

The morphologies of galaxies formed in the simulations with 
shock-induced bursts are rather sensitive to the value of the burst 
threshold A'burst, as shown in §4.4, because the allowed range for 
this parameter is very narrow. With too low a threshold, the star for¬ 
mation is dominated by the burst mode at all times and vice versa. 
As a result, the allowed range for Aburst is less than an order of 
magnitude in our simulations. Within this range, the behaviour of 
the shock-burst models is easily understood. Lower threshold val¬ 
ues suppress star formation more effectively at high redshift. In¬ 
terestingly, the three shock-burst models have almost the same star 
formation rate at low redshift, when no major mergers occur. Since 
our current modelling of shock-induced bursts is not resolution in¬ 
dependent as we discussed in §2.2.4, we will explore alternative 
formulations in future work. 

Several important tests of our model are possible and we in¬ 
tend to pursue these in future work. For example, we have followed 
in detail the evolution of the metallicity of gas and stars in our simu¬ 
lations, keeping track of metals produced in the quiescent and burst 
modes. Comparing the resulting distributions of heavy elements in 
the different dynamical components of the simulated galaxies with 
observational data will be an important test of the validity of our 
feedback models. A key question is how often do shock-induced 
bursts like those we have modelled lead to the formation of disc- 
dominated galaxies. Answering this question will require simulat- 
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ing galaxy formation in large volumes. The results presented in this 
paper encourage us to believe that such demanding calculations are 
worth pursuing. 
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